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INTRODUCTION 


ATMOSPHERIC  RADIATION 

The  funds  for  the  atmospheric  radiation  portion  of  the 
contract  were  exhausted  as  of  the  previous  report;  however, 
work  has  continued  with  other  support.  A  summary  of  the  prog¬ 
ress  made  follows. 

Major  effort  was  directed  towards  an  extensive  study 
of  solar  radiation  fluxes  in  the  presence  of  typical  summer 
Arctic  stratus  clouds.  A  paper  on  this  subject  was  delivered 
by  Dr.  Warren  J.  Wiscombe  at  the  24th  Alaskan  Science  Confer¬ 
ence  in  Fairbanks,  Alaska  (August  15-17,  1973).  A  consider¬ 
ably  extended  and  generalized  version  of  this  paper  has  been 
submitted  for  publication  in  the  proceedings  of  that  confer¬ 
ence. 

Massive  code  modifications  were  made  in  ATRAD  (the 
atmospheric  radiation  computer  model)  which; 

(1)  introduced  a  specular  reflection  option 
(important  for  calculations  involving 
smooth  or  slightly  roughened  sea  sur¬ 
faces)  ; 

(2)  allowed  the  use  of  a  grey-body  top 
boundary  with  arbitrary  temperature  and 
emissivity,  so  that,  for  example,  long¬ 
wave  (IR)  calculations  need  not  include 
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the  whole  atmosphere  but  may  be  termi¬ 
nated  at  a  cloud  base; 

(3)  streamlined  the  basic  Grant-Hunt  al¬ 
gorithm  computation  significantly; 

(4)  allowed  the  code  to  cycle  over  any  one 
of  the  following  parameters  —  sun 
angle,  surface  reflectivity,  surface 
temperature,  top-boundary  temperature. 

The  cycling  feature  in  (4)  is  a  tremendous  aid  in  making  ex¬ 
tensive  parameter  studies  such  as  were  required  for  the  Arctic 
stratus  problem,  and  was  achieved  with  virtually  no  increase 
in  computing  time  over  previous  calculations  involving  a 
single  parameter.  This  was  possible  because  most  of  the 
calculations  in  ATRAD  do  not  depend  on  a  particular  parameter, 
and  thus  it  was  possible  to  cycle  only  those  pares  of  the 
calculation  depending  on  the  parameter  in  question.  Since  any 
one  of  four  parameters  may  be  cycled,  the  code  modifications 
were  necessarily  complex  and  time-consuming. 

Finally,  a  seminar  "Comparisons  of  a  Detailed  Radia¬ 
tion  Model  with  the  Radiation  Subroutine  of  the  Mintz-Arakawa 
General  Circulation  Model"  was  prepared  and  presented  as  part 
of  the  weekly  seminar  series  of  the  UCLA  Meteorology  Depart¬ 
ment. 

OROGRAPHIC  EFFECTS  ON  GLOBAL  CLIMATE 

The  results  reported  below  pertain  to  the  continuation 
of  numerical  investigations  of  meso-scale  momentum  transfer 
in  the  atmosphere.  Investigations  can  be  separated  into  two 
parts:  the  transient  and  the  steady-state. 
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The  steady-state  part  has  culminated  in  the  develop¬ 
ment  of  two-  and  three-spatial  dimensional  codes  based  on  the 
formulation  of  Bretherton.  ^  Additional  analysis  of  the 
Bretherton  equations  is  found  in  Section  2,  and  a  3-D  study  of 
the  wave  drag  over  the  Sierra  Nevada  range  of  California  is 
described  in  Section  3.  The  steady-state  computer  code  is 
divided  into  three  parts  to  more  efficiently  conduct  param¬ 
eterization  studies.  This  work  is  described  in  Section  4. 

Transient  phenomena  have  been  investigated  with  the 
2-D  code  HAIFA,  and  the  3-D  code  STUFF.  The  latter  has  been 
optimized  in  terms  of  speed  of  execution  and  of  fast  storage 
requirements.  A  complete  description  of  the  code  optimization, 
along  with  a  discussion  of  the  numerical  techniques  employed 
by  STUFF,  is  given  in  Part  II  of  this  report. 

During  the  next  contract  period  major  emphasis  «il3  be 
placed  on  application  of  the  transient  and  steady-state  codes 
to  the  development  of  heuristic  expressions  characterizing 
meso-scale  momentum  transfer  processes  induced  by  global 
topography. 
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1.  INTRODUCTION 


/ 

By  the  close  of  the  last  contract  period,  the  devel¬ 
opment  of  a  3-D,  linear,  steady-state  code  based  on  the  Brether- 
ton ^  formulation  had  been  completed.  An  outline  of  the  code 
was  also  presented.  Since  then,  several  subtle  points  of 
Bretherton's  paper  have  been  analyzed.  A  discussion  of  these 
is  presented  in  Section  2.  Section  3  presents  the  results  of 
a  3-D  steady-state  study  of  the  Sierra  Nevada-Owen' s  Valley 
region  of  north-central  California  where  lee-wave  measurements 
have  been  carried  out. 

Analysis  of  the  3-D  code  during  this  contract  period 
has  shown  that  it  might  be  significantly  optimized  by  reorga¬ 
nization  and  recoding.  The  results  of  this  task  are  presented 
in  Section  4.  A  discussion  of  the  numerical  method  is  in¬ 
cluded.  Finally,  a  discussion  of  plans  for  future  work  is 
presented  in  Section  5. 
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2.  ADDITIONAL  ANALYSIS  OF  THE  BRETHERTON  FORMULATION 

The  Eq.  (54)  of  Bretherton/^  on  which  the  3-D  steady- 
state  evaluation  of  orographic  gravity  wave  drag  is  based,  was 
presented  in  his  original  paper  substantially  without  deriva¬ 
tion.  Subsequently,  several  of  the  steps  in  the  derivation 

[21 

were  presented  in  a  recent  progress  report  on  calculations 
usin  the  Bretherton  formulation.  Several  points  in  the  Bre- 
therton  paper,  however,  which  were  not  examined  at  that  time, 
subsequently  have  been  found  to  contain  subtitles.  In  addi¬ 
tion,  errors  in  the  text  have  been  recognized  and  corrected. 
Consequently,  we  briefly  record  below  the  salient  arguments 
required  to  complete  the  Bretherton  formulation. 

2.1  CHOICE  OF  VERTICAL  VELOCITY  AT  UPPER  BOUNDARY 

The  solution  cf  the  vertical  velocity  equation  re¬ 
quires  boundary  conditions  to  be  imposed  both  at  the  ground 
and  at  the  top  of  the  atmosphere.  In  order  to  perform  a 
marching  calculation,  a  specific  value  of  the  complex  vertical 
velocity  is  assigned  at  the  top  of  the  atmosphere.  The  result¬ 
ing  wave  drag  is  independent  of  this  value,  a  result  which  we 
wish  to  demonstrate  below.  In  order  to  do  so,  we  evaluate  the 
quantity  F  containing  all  of  the  vertical  velocity  terms. 
According  to  Bretherton,  we  must  consider  the  two  cases  for 
which  different  boundary  conditions  at  the  top  of  the  atmos¬ 
phere  z=H  are  prescribed. 
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Before  considering  these  two  cases,  however,  we 
develop  a  general  solution  of  the  Scorer  equation  with  which 
we  can  more  easily  consider  the  boundary  conditions.  The  real 
and  imaginary  parts  of  the  vertical  velocity  each  obey  the 
same  linear  second  order  equation 

fjT  +  U2-<2)u  =  0  .  (1) 

A  general  sciut-.on  of  this  equation  contains  two  coef¬ 
ficients 


u 


a-^  +  a2u2 


(2) 


where  u^^  and  u2  are  linearly  independent.  We  choose  the 
functions  u1  and  u2  to  be  solutions  of  Eq.  (1)  having  the 
following  boundary  conditions: 


dM 

du2\  - 1 
**-)  h'  1 


(3) 


As  mentioned  above,  the  real  and  imaginary  parts  of  w  ,  the 
vertical  component  of  the  perturbation  velocity,  each  have  the 
form  of  Eq .  (2)  . 

We  now  discuss  the  simpler  of  the  two  cases  in  which 
k2  >  S,2(h)  ,  for  which  we  have  trapped  wave  solutions  and  the 
contributions  to  the  Reynold's  stress  are  discrete  (as  given  by 
F.q.  (61)  of  Bretherton)  .  The  boundary  condition  at  the  top  of 
the  atmosphere  z=H  in  this  case,  is  Eq.  (50a)  of  Bretherton: 
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dw 

dl 


=  -  /fT 


i2  (H)ti 


Assuming  an  initial  value  of  w  at  z=H  of  ft(H)  = 
uQ  +  ivQ  ,  the  initial  values  of  the  derivatives  are 

ar)H  =  -  - 

sr)H  -  -  •/^rmv0  . 

We  can  now  determine  the  coefficients  a1  and  a2  in 

Eq.  (2): 


WR  =  Uo(ul  "  Sk2-12  (H)u2) 
WI  =  vo(ul  “  /<2-^2  (H)u2) 


(4) 


Thus,  we  find  that  the  real  part  of  the  velocity  wD  and  the 
imaginary  part  w^  are  proportional  to  each  other  at  every 
altitude  z  .  Consequently,  it  is  only  necessary  to  perform 
one  integration  of  the  Scorer  equation  to  evaluate  either  one 
of  them. 


It  is  now  possible  to  construct  the  term  containing 
all  of  the  velocity  dependence  of  the  Reynold's  stress,  given 
by  Bretherton's  Eq.  (60); 
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In  terms  of  the  expressions  of  Eq.  (4),  we  obtain 


/ 


F (k) d< 


TT 

2k 


y<2-i2  (h) 


z=0 


-  A2-l2  (H)u2)  2  dz 


(5) 


The  salient  feature  of  Eq.  (5)  is  that  the  valuer,  cf 
UQ  and  vQ  do  not  appear;  the  trapped  wave  contribution  to 
the  Reynold's  stress  is  independent  of  the  assumed  velocity 
amplitude  at  z=H  . 

For  the  case  of  waves  which  leak  into  the  stratosphere, 
coriesponding  to  k2  <  5, 2  (H)  ,  the  boundary  condition  (given 
by  Eq.  (50b)  of  Bretherton)  is 

§§  =  +  i A2  (H)  -k2  sgn(Un)w  . 

Substituting  the  assumed  boundary  values  at  z=H  the  deriva¬ 
tive  condition  becomes 

=  -  A2  ( H /  —  k 2  sgn(Un)vQ  , 

=  A2  (H)  -<2  sgn(Un)uo 


The  general  solutions  satisfying  these  boundary  conditions  are 

*r  =  uiuo  "  'H)-k2  sgn(Un)u2vo  , 

/ -  (6) 

Wj  =  u1vo  +  /£2(H)-k2  sgn(Un)u2uo  . 
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In  contrast  to  the  first  case,  these  two  solutions  are  linearly 
independent  and  require  that  two  integrations  of  the  Scorer 
Equation  be  performed. 

In  order  to  evaluate  the  term  containing  all  of  the 
dependence  on  the  vertical  velocity  in  the  Reynold's  stress, 
we  compute  the  quantity  F  (see  Bretherton's  Eq.  (52)): 


F 


( dw 


dWA* 

dz 


6«nft*(0) 


(7) 


where  w*  denotes  the  complex  conjugate  of  w  .  The  numerator 
of  F  is  independent  of  altitude  and  can  be  evaluated  at  z=H 
from  the  boundary  values  given  above.  The  result  is 


=  A2  (H)  -k2  sgn(Un)  fo(H)  0*  (H) 
=  A2  (H)  -k2  sgn  (Un)  (u2+v2 ) 


Using  this  result,  and  forming  the  denominator  of  Eq.  (7)  from 
the  solutions  of  Eq.  (6)  evaluated  at  z=0  ,  we  obtain  for  F 


A2  (H)  -k 2  sgn  (U  ) 

F  -  - — - 

u2  (0)  +  (S,2  (H)-k2)u2 (0) 


(8) 


This  expression  is  also  independent  of  the  assumed  boundary 
values  uQ  and  vq  .  Consequently,  we  have  shown  that  the 
Reynold's  stress  does  not  depend  cn  the  chosen  values  of 
vertical  velocity  at  the  top  of  the  atmosphere. 


10 


SSS-R-74-2023 


2 . 2  REDUCTION  OF  DOMAIN  OF  INTEGRATION 

This  result/  which  reduces  the  amount  of  calculation 
in  forming  the  Reynold's  stress  components  by  a  factor  2,  is 
incorrectly  justified  by  Bretherton.  In  fact,  as  will  be 
shown  below,  the  simplification  results  from  the  integrand 
being  the  same  at  the  wave  numbers  (k,£)  and  (-k,-£)  rep¬ 
resenting  a  reflection  in  the  origin  of  wave  number  space, 
rather  than  as  stated  by  Bretherton  that  (k,4>)  is  the  same 
as  (k,-<|>)  . 

Several  factors  enter  in  the  Reynold's  stress  inte¬ 
grand;  we  examine  the  behavior  of  each  of  these  when  -uhe 
transformation  (-k,-*)  ■+  (k,fc)  is  carried  out.  First,  we 
examine  the  differential  equation  for  the  vertical  velocity 
and  its  boundary  conditions. 


Kz(-k,-£)  =  k2+l2  =  k2  (k,  i.) 


u n(-k,-£)  =  ^(Uk+V£)  =  -  U  (k,  £) 

,  N2  d2U 

a  (-k,-£)  =  jjj-  -  — £2  =  £*(k,£) 

n  n 


(9) 


i.e.,  the  Scorer  parameter,  being  an  even  function  of  IJ 

n  ' 

is  invarient  under  the  transformation.  Consequently,  the 
Scorer  equation  and  its  elementary  solutions,  depending  only 
on  fi.2-<2  ,  are  also  invarient.  Using  this  result,  which 
establishes  that  u^  and  U2  of  Eq.  (3)  are  invarient,  and 
considering  the  transformation  properties  of  Eq.  (8)  for  F 
we  find  that 


F(-k,-J>)  =  -  F(k,Jl) 


(10) 
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We  now  consider  the  topographic  factor  of  the  Reynold' 
stress  given  by  Eq.  (47)  of  Bretherton. 


fi(-k,-U)  =  ^  £  £  h(x,y>  e-U-kx-ly)  gxdy 

=  fi* (k# 1)  . 

Forming  the  spectrum  function  A  =  fi  fi*  ,  we  obtain  the 
transformation  property 

A(-k,-£)  =  ~~  h  ( — k ,  —  £ )  h*  (-k,-£)  =  A(k,Jt)  . 


These  quantities  can  be  combined  to  form  the  integrands  of  the 
two  horizontal  components  of  the  Reynold's  stress: 


puw' 

pvw. 


P  U2(0)  k2Af| 


o  n 


!cos4>  ■ 
sin<f) , 


Po0*<0> 


kAF 


The  transformation  (-k,-M  +  (k,A)  is  seen  to  leave  the  inte 
grand  invarient: 


puw(-k,-£)  =  puw(k,£) 
pvw (-k,-£)  =  pvw (k, Z) 

Consequently,  the  result  given  in  Bretherton' s  Eq.  (54)  is 
confirmed;  namely  that 

/*2tt  r  tt 

I  d<p  =  2  i  dtp  . 

Jo  Jo 
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3.  COMPARISON  OF  TRANSIENT  AND  STEADY-STATE  WAVE  DRAG 

3.1  INTRODUCTION 

At  the  close  of  the  last  contract  period,  two  investi¬ 
gations  of  wave  drag  phenomena  associated  with  the  orography 
of  the  Sierra  Nevada  had  been  completed.  These  were  the  2-D 
transient  study  presented  in  the  preceding  report,121  and  the 
2-D  linear  steady-state  study  discussed  in  another  report. ^1 
In  the  transient  study,  atmospheric  flow  over  the  Sierra 
Nevadas  was  simulated,  resulting  in  the  time— dependence  of 
wave  drag  values.  Additionally,  the  effects  of  moisture  on 
wave  drag  were  evaluated.  A  steady-state  calculation  was  also 
performed,  from  which  a  small-scale  parameterization  using 
some  of  the  important  dynamic  quantities  was  developed.  The 
results  of  this  calculation,  where  applicable,  were  compared 
to  the  transient  results. 

3.2  THREE-DIMENSIONAL  STEADY-STATE  STUDY 

After  the  completion  of  the  3-D  steady-state  code,  two 
test  problems  were  run.  The  first  investigated  a  generaliza- 
tion  of  the  triangular  mountain  of  the  previous  report^1  by 
calculating  a  triangular  ridge  128  km  in  length.  The  second 
calculation  was  applied  to  the  accual  Sierra  Nevada—Owen1 s 
Valley  topography  described  below. 

The  Sierra  Nevada  is  a  single  unbroken  range  approxi¬ 
mately  400  miles  in  length  and  50  to  80  miles  wide.  At  the  • 
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Owen's  Valley  we  find  the  mam  crest  running  nearly  North- 
South.  Elevations  along  the  crest  are  around  12,000  ft.  with 
numerous  peaks  rising  above  14,000  ft.  The  eastern  scarp  is 
abrupt  and  quite  straight.  Owen's  Valley  is  to  the  east  at  an 
average  elevation  of  4,000  ft.  and  a  nearly  uniform  width 
of  approximately  15  miles.  The  east  wall  of  Owen's  Valley  is 
a  fault  block  range  called  the  Inyo  Mountains  to  the  south 
and  the  White  Mountains  to  the  north.  The  elevation  of  the 
Inyo  Mountains  is  7,000-11,000  ft.  while  the  White  Mountains 
rival  the  Sierras  in  height. 

Figures  1  and  2  show  the  topography  of  the  Owen's 
Valley  region.  Figure  1  is  a  contour  map  of  the  data  used  in 
the  second  3-D  calculation.  Figure  2  shows  three  east-west 
transects  at  the  indicated  latitudes.  The  asterisk  indicates 
the  position  of  Owen's  Valley.  Both  figures  were  made  from 

data  derived  from  the  5'  x5'  topography  tapes  described  in 
Section 


The  steady-state  runs  utilized  the  wind  and  tempera¬ 
ture  profile  data  obtained  from  the  Merced  weather  station  and 
employed  in  the  previous  studies .  [2 ' 3 ]  The  calculations  indi¬ 
cate  a  principal  trapped  wave  having  a  wavelength  of  approxi¬ 
mately  19  km.  This  agrees  well  with  the  value  observed  in  the 
Owen's  Valley  under  the  same  meteorological  conditions  of 
18-20  km.  The  wave  drag  results  are  shown  superimposed  on  the 
transient  HAIFA  results  in  Figure  3.  Two-dimensional  steady- 
state  calculation  results  utilizing  the  east-west  transects  of 
Figure  2  are  also  presented.  As  can  be  seen,  the  results  of 
the  3-D  runs  agree  quite  closely,  revealing  a  possible  insen¬ 
sitivity  to  the  topography. 

Ihe  2-D  runs,  however,  differ  by  as  much  as  a  factor 
of  3.  Since  the  real  topography  calculations  and  idealized 
topography  3-D  calculations  agree  closely,  one  may  assume  that 
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Figure  2 


X  |k«) 


2-D  cross-sections  of  Sierra  topography  from  120°W 
through  116°50'W  longitude.  Asterisk  indicates 
Owen's  Valley. 
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the  idealized  topography  is  a  good  characterization  of  the 
real  topography.  The  differences  between  the  2-D  and  3-D  re¬ 
sults  can  thus  be  attributed  to  the  sensitivity  of  Reynold's 
stress  to  the  topographic  variability  exhibited  in  Figure  2. 

These  results  provide  some  indirect  evidence  of  the 
accuracy  with  which  topographic  data  must  be  represented.  The 
question  of  sensitivity  to  the  representation  of  topographic 
data  will  be  examined  in  greater  detail  during  the  next  con¬ 
tract  period. 

It  is  of  interest  to  examine  the  wavenumber  dependence 
of  the  atmospheric  factor  F  (the  spectral  Reynold's  stress 
corresponding  to  unit  topographic  spectrum  function) .  This 
quantity  can  indicate  to  us  what  spectral  regions  are  most 
significant  in  serins  of  both  the  continuous  and  trapped  wave 
Reynold's  stress  contributions.  We  examine  below  this  quan¬ 
tity  for  the  atmospheric  parameters  of  the  Merced  problem. 

The  distribution  of  the  trapped  waves  is  seen  in  Figure  4. 

The  locus  for  a  particular  trapped  wave  is  generated  by  exam¬ 
ining  the  solution  of  the  Scorer  equation  for  each  of  the 
plotted  points.  The  locus  is  then  constructed  by  connecting 
points  whose  solutions  have  identical  numbers  of  nodes.  Cur¬ 
rent  A(f>,A<  values  allow  the  identification  of  three  trapped 
waves.  These  waves  are  modest  contributors  to  the  drag  for 
the  Sierra  Nevada  problem.  Apparently,  other  waves  which  are 
not  adequately  resolved  in  the  figure  can  be  substantial  con¬ 
tributors  as  well.  Some  of  the  individual  points  not  associ¬ 
ated  with  an  established  locus  have  been  found  to  contribute 
more  than  those  found  on  the  loci.  It  appears  as  though 
integrations  having  higher  resolution  in  Ak,A<{>  are  required 
to  determine  these  trapped  wave  contributions  more  fully. 

One  worrisome  feature  of  the  analytic  expression  for 
the  trapped  wave  contribution  to  the  drag  is  that  waves  whose 
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Figure  4  Trapped  wave  distribution  in  <r<t>  space.  Points 
not  falling  on  curves  represent  contributors  to 
unresolved  trapped  waves.  N  is  the  number  of  nodes 
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loci  are  essentially  radial  in  <,<p  space  may  be  inaccurately 
integrated  in  the  wave  drag  calculation.  This  is  due  to  the 
fact  that  an  increment  of  A<{>  may  miss  the  trapped  wave  par¬ 
tially  or  even  completely.  Some  of  the  loci  of  the  waves  of 
Figure  4  exhibit  quite  nearly  radial  behavior  and  consequently 
some  of  the  contribution  may  be  lost.  A  change  to  integrating 
over  k,i  instead  of  <  <£  as  independent  variables  might 
alleviate  this  since  the  spatial  resolution  would  be  more 
easily  controlled  at  greater  distances  from  the  origin. 

Figure  5  is  a  contour  plot  of  the  continuous  contribu¬ 
tion  to  the  F  distribution  for  the  Sierra  Nevada  problem. 
Both  the  conti nuous  and  the  trapped  waves  exhibit  symmetry 
about  <f>=0°  in  this  case.  The  reason  for  this  can  be  seen  in 
the  symmetry  'if  the  Scorer  parameter  and  the  boundary  condi¬ 
tion,  which  in  turn  implies  symmetry  in  the  solutions  of  the 
Scorer  equation.  For  the  Sierra  problem  Un  is  given  by 

Un  =  U  ( z)  cos<(> 

Since  the  symmetry  of  the  Scorer  parameter  depends  on  the  sym¬ 
metry  properties  of  Un  ,  we  see  that  the  Scorer  parameter 
will  be  symmetric  about  <f>=0°  . 

The  F  values  in  Figure  5  remain  large  even  for  < 
values  quite  close  to  Jt(H)  •  In  fact,  for  some  4>  ,  the  F 
values  actually  increase  until  <  =  l (H)  .  This  behavior  can 
be  explained  by  examining  the  expression  for  F  : 

F  =  2 /i2  (H)  -k2  sgn(Un)/w*  (0)w(0)  . 

One  sees  that  the  K-dependence  of  F  is  predominantly  con¬ 
trolled  by  the  denominator  for  <  <<  1(H)  .  For  values  of  k 
near  S,(H)  ,  however,  the  behavior  of  F  can  be  affected  by 
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both  numerator  and  denominator.  If  the  denominator  is  not 
singular  near  k=£(H)  7  will  decay  monotonically  to  zero  as 

k  approaches  £(H)  .  However,  if  the  denominator  has  a  zero 
near  l (H)  ,  a  resonant  behavior  can  occur  in  F. 

This  behavior  corresponds  to  the  existence  of  a  trap¬ 
ped  wave  near  MH)  .  While  this  situation  was  not  expected 
to  occur  very  often,  the  Sierra  Nevada  calculations  reveal 
several  such  cases  as  shown  in  Figure  6.  The  methods  used  in 
the  numerical  integration  currently  do  not  take  account  of 
this  occurrence  accurately;  an  improved  integration  routine 
will  be  developed  to  deal  with  this  case. 

Apparently,  this  near  coincidence  of  a  resonance  with 
the  locus  of  £ (H)  did  not  take  place  in  the  example  calcu¬ 
lated  by  Bretherton.  As  a  consequence,  he  did  not  find  it 
necessary  to  take  special  care  with  the  integration  routine 
in  this  region.  As  additional  examples  are  investigated  it 
is  to  be  expected  that  other  special  cases  will  be  found 
which  also  will  require  special  consideration. 
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Figure  6  F(k)  distribution  for  various  <p  angles. 
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4.  FURTHER  CODE  DEVELOPMENTS 


4.1  INTRODUCTION 

When  the  steady-state  codes  became  operational,  various 
test  problems  were  run  to  determine  their  response  to  differ¬ 
ent  problem  configurations.  The  comparison  of  the  results  of 

f  21 

a  3-D  run  and  the  previously  reported  1  HAIFA  study  of  the 
Sierra  Nevada  problem  were  presented  in  Section  3.  They  indi¬ 
cate  that  similar  Reynold  stress  values  are  obtained  by- these 
two  markedly  different  codes. 

On  completion  of  this  comparison,  work  was  begun  to 
improve  the  efficiency  of  the  steady-state  code,  preparatory 
to  a  large  number  of  survey  calculations.  This  program  is 
discussed  more  fully  in  Section  5.  In  summary,  the  code  is 
segmented  such  that  the  spectrum  function,  F  factor  and  stress 
integrals  are  calculated  independently  of  each  other,  but  ex¬ 
change  data  through  files  established  on  mass  storage  devices. 
This  procedure  makes  efficient  use  of  fast  core  storage  and 
decreases  computer  costs.  At  the  same  time,  much  greater  ef¬ 
ficiency  is  achieved  in  performing  a  matrix  of  calculations. 
Since  this  modulization  has  resulted  in  major  design  changes 
in  the  codes,  a  description  of  the  new  routines  follows.  The 
topography  data  files  have  also  undergone  major  modification 
and  streamlining.  These  are  also  discussed  below. 
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4.2  TOPOGRAPHY 

Topography  data  were  obtained  from  the  Defense  Mapping 
Agency  (DMA)  and  are  contained  on  two  sets  of  magnetic  tapes. 
The  first  set  contains  elevations  area-averaged  over  30'x30' 
rectangles  for  the  entire  globe.  The  second  set  contains  ele¬ 
vations  area-averaged  over  5' *5'  rectangles  for  most  of  North 
America  and  Europe. 

An  examination  of  these  tapes  indicated  that  a  re¬ 
structuring  of  the  data  format  would  allow  much  faster  access 
to  particular  pieces  of  information.  As  a  result,  the  reading 
of  all  available  30'  data  now  requires  30  seconds  of  CPU  time 
as  compared  to  600  seconds  previously. 

The  new  tapes,  resulting  from  the  above  restructuring, 
are  arranged  in  r  common  format  for  both  the  30'  and  5'  data 
sets.  The  30'  set  requires  two  tapes,  one  for  the  northern 
hemisphere  and  one  for  the  southern  hemisphere.  The  5'  set 
consists  of  one  tape  covering  the  northern  hemisphere.  The 
data  are  written  in  logical  records  having  1444  word  length 
for  30'  data  and  4324  words  for  5'  data.  Each  logical  record 
represents  topographical  data  at  a  particular  latitude  for  a 
full  360°  of  longitude.  These  records  are  written  on  the 
tapes  sequentially  starting  at  the  northernmost  latitude  of 
each  hemisphere.  At  the  head  of  each  logical  record  are  writ¬ 
ten  the  starting  latitude  and  longitude  of  the  strip,  where 
longitude  is  measured  east  of  Greenwich.  Associated  with  each 
data  point  is  a  code  number  which  contains  information  about 
the  nature  of  the  topography;  i.e.,  whether  sea  bed,  lake,  all 
land,  etc.  This  coded  information  is  the  same  as  that  in  the 
DMA  description  of  their  tapes.  In  addition  to  the  code  data, 
however,  another  indicator  has  been  added  to  denote  missing 
data.  It  was  found  that  occasional  gaps  in  data  occurred  on 
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the  DMA  tapes.  Since  the  Fourier  transform  routine  requires 
complete  data,  a  flag  is  set  whenever  the  missing  topography 
is  encountered  and  the  run  is  terminated.  Usually,  this  re¬ 
quires  obtaining  the  data  from  other  sources  and  inserting  it 

into  the  data  tape.  The  case  can  then  be  rerun  in  the  usual 
manner. 

4.3  THE  SPECTRUM  FUNCTION  CODE 

4.3.1  The  Equations 

The  spectrum  function  code  performs  the  primary  task 
of  determining  a  spectrum  function  distribution  in  wavenumber 
space  associated  with  a  Fourier  representation  of  the  hori¬ 
zontal  dependence  of  the  topography.  Bretherton defines  a 
spectrum  function,  A; 

A (k,  1)  =  ill  fi*  fi  (] 

where  X  and  Y  represent  the  spatial  extent  of  the  topo¬ 
graphic  region  in  the  east-west  and  north-south  directions, 
is  the  Fourier  transform  of  the  surface  height  and  the  asterisk 
indicates  conjugation.  The  Fourier  transform  is  defined  as 


h  (k ,  1) 


1  /*y  -i(kx+ly) 

47F  Jo  J  h(x,y)e  dxdy 


(2) 


These  two  equations  are  solved  in  Module  1. 

4.3.2  Numerical  Method 

The  approach  taken  in  Module  1  is  as  follows:  First, 
surface  height  data  are  obtained  from  the  appropriate  data 
tape.  Next,  these  data  are  Fourier  transformed  to  obtain  the 
spectrum.  Finally,  the  spectrum  function  distribution  is  cal¬ 
culated  according  to  Eq.  (1). 
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The  topography  data  are  composed  of  m*n  data  "points" 
falling  within  the  rectangle  specified  by  the  user  at  problem 
generation  time.  It  must  be  remembered  that  each  "point"  is, 
in  fact,  an  average  value  to  be  associated  with  a  small  rec¬ 
tangle  about  the  point  in  question.  Thus,  the  surface  height 
may  be  thought  of  as  a  step-function  representation  of  the 
topography.  A  sample  5*4  grid  shown  in  Figure  7  demonstrates 
this  more  clearly.  Each  of  the  numbers  appearing  within  the 
cells  represents  the  average  height  within  the  cell. 


[2 


The  nature  of  the  data  representation  has  implications 
for  the  Fourier  transform.  First,  the  discrete  data  can  be 
represented  by  a  discrete  number  of  spectral  components. 

Second,  the  resolution  in  spectral  space  is  a  function  of  the 
resolution  in  real  space.  As  mentioned  in  a  previous  report, 
this  lack  of  spectral  resolution  can  affect  the  wave  drag  cal¬ 
culation  through  rendering  A  uncertain.  Particularly  at 
trapped  wavenumbers  this  uncertainty  results  in  a  corresponding 
uncertainty  in  the  trapped  wave  drag  contribution.  The  use  of 
5'  data  allows  resolution  of  wavelengths  of  the  order  of  10  km, 
which  are  comparable  with  the  wavelengths  of  the  predominant 
trapped  waves.  We  propose  to  examine  the  desirability  of  using 
even  more  highly  resolved  data. 


The  fact  that  the  topography  can  be  considered  discrete 
is  fortunate  from  a  calculational  point  of  view,  however.  Cal¬ 
culation  of  the  Fourier  transform  requires  a  large  amount  of 
time  due  to  the  large  number  of  trigonometric  evaluations  and 
manipulations  of  the  data.  The  fast  Fourier  transform  (FFT) , 
which  is  an  algorithm  to  optimize  the  direct  calculation  of  the 
Fourier  transform,  is  used  to  calculate  the  topography  spectrum 
function.  The  algorithm  requires  that  the  function  be  defined 
at  a  discrete  set  of  points  spaced  at  equal  intervals  in  each 
dimension.  A  FORTRAN  subroutine  based  on  the  Cooley  Tukey 
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algorithm  was  obtained  from  the  University  of  California  at 
San  Diego  and  has  been  used  satisfactorily  in  the  applications 
described  below. 

The  finite  Fourier  transform  requires  some  adaptation 
before  it  is  applicable  to  the  topography  data.  This  can  be 
seen  from  the  finite  difference  expression  for  the  Fourier 
integral  which  is  given  by 


m-1  n-1 


=  jp  ^  h  (sAx,uAy)  e 


-i (ksAx+luAy) 


AxAy 


s=0  u=0 


(3) 


The  FFT  routine  calculates 


m-1  n-1  _  .  /rs  tu\ 

; .  =  i-  y  y  .  e‘2ni^  +  “) 

rt  mn  /  j  /  j  su 

s=o  u=0 


(4) 


We  see  then  that  if  we  identify  a  with  the  height  function 

o  U 

as  follows 


a 


su 


mAx  nAy 
4tt2 


h (sAx,uAy) 


(5) 


and  define 


2tTr 

mAx 


0  <  r 


(6) 


2-tt  (r-m) 
mAx 


<  m 
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1  = 


2jrt 

nAy 


0  1  t  <  £ 


(7) 


_  2 tt  (t-n) 
nAy 


n 

2 


<  t  <  n 


then  the  FFT  routine  will  return  the  Fourier  transform  as  de¬ 
fined  in  Eq.  (2)  .  The  wavenumbers  associated  with  the  trans¬ 
form  are  given  by  Eqs.  (6)  and  (7). 

The  DMA  data  tapes,  however,  do  not  contain  values 
associated  with  boundary  points  of  the  area  in  question.  In¬ 
stead,  the  data  are  associated  with  averages  over  elementary 
rectangles  ard  lie  in  the  interior  of  the  area  as  illustrated 
by  Figure  7.  It  is  also  necessary  that  the  number  of  data 
points  along  each  spatial  dimension  be  an  integer  power  of  2. 

In  order  to  take  account  of  this  displacement  of  the  data  from 
the  boundary  points,  a  redefinition  of  the  Cooley-Tukey  coeffi¬ 
cients  is  required: 


asu  *  — 4^y  ht(s-H,)Ax,  (u+Hliy] 


(8) 


Following  the  calculation  of  the  FFT  coefficients  a  subsequent 
calculation  to  form  the  topography  transform  is  carried  out: 

A  „  -h i  (kAx+lAy) 

hk£  =  art  e  (9) 

This  transformation  of  terms  allows  one  to  approximate  the 
integral  of  Eq.  (2)  without  having  to  use  boundary  points. 

The  steps  of  the  spectrum  function  calculation,  conse¬ 
quently,  are  as  follows:  The  topography  data  are  interpolated 
to  obtain  new  values  corresponding  to  the  next  higher  integer 
power  of  2  equally  spaced  points  in  each  direction.  These 


30 


SSS-R-74-2023 


values  are  pre-processed  according  to  Eq.  (8) .  The  FFT  is 
executed  and  the  resulting  array  is  then  processed  as  in  Eq. 

(9)  to  obtain  the  Fourier  transform  of  the  height.  The  spec¬ 
trum  function  is  then  straightforwardly  calculated  according 
to  Eq.  (1/ . 

The  spectrum  function,  together  with  the  associated 
k,l  wavenumbers,  is  written  into  a  file  for  future  use.  The 
final  task  of  Module  1  is  to  print  a  tabulation  of  A(k,l) 
Contour  maps  of  the  topography  and  spectrum  function  are  also 
plotted.  Additionally,  a  contour  map  of  the  log  of  the  spec¬ 
trum  function  is  made. 

4.4  SPECTRUM  FUNCTION  CALCULATIONS 

Since  the  spectrum  function  code  was  completed  late  in 
the  last  contract  period,  only  a  few  example  calculations  of 
the  topography  function  A  have  been  completed.  The  study 
consisted  of  an  examination  of  ten  geographic  areas  in  the 
United  States  differing  widely  in  topography  and  location. 
Topography  data  from  the  DMA  magnetic  tapes  having  5'  resolution 
were  used  and  the  resulting  spectrum  functions  were  calculated. 
These  same  topographic  data  were  then  resolved  to  10',  by  com¬ 
bining  5'  values,  and  the  spectrum  functions  were  determined 
again.  The  resulting  set  of  data  can  be  investigated  to  deter¬ 
mine  the  behavior  of  A  as  a  function  of  differing  topographies 

and  also  determine  the  sensitivity  of  A  to  topographical  reso¬ 
lution. 

A  representative  sample  of  terrains  was  selected;  these 
data  consisted  of  contiguous  rectangles  within  a  strip  of  land 
across  the  U.S.  from  34°20’  to  37°0'  north  latitude  and  from 
94 ° 20 1  to  121°  west  longitude.  This  strip  includes  the  rugged 
topography  of  the  High  Sierra  and  Rocky  Mountains,  the  flat 
plains  of  the  midwest,  and  the  Appalachian  mountains  of  the  east. 
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The  strip  was  composed  of  10  contiguous  topographical  grids 
2°40 '  on  a  side  (corresponding  to  32x32  data  points  for  each 
grid  at  5  resolution) .  This  set  of  grids  constitutes  a  pre¬ 
liminary  cross-section  of  characteristic  topographies. 

Since  the  study  is  still  in  progress,  only  a  brief 
examination  of  a  sample  grid  is  currently  available.  Figures 
8-11  show  contour  maps  of  one  geographic  area  and  its  associ¬ 
ated  spectrum  function  located  in  the  rectangle  118° 20'  to 
121°  west  longitude,  34°20'  to  37°0'  latitude  at  both  5'  and 
10'  resolution.  This  area  corresponds  to  the  first  grid  in 
the  set  of  grids  and  encompasses  portions  of  the  San  Joaquin 
Valley ,  Sierra— Nevada  mountains  and  the  coastal  ranges  of 
California.  Figure  8  presents  the  contour  map  at  5'  resolu¬ 
tion  and  Figure  9  is  the  corresponding  calculated  spectrum 
function..  Figure  10  presents  the  contour  map  for  the  10' 
resolution  data.  Figure  11  is  the  corresponding  spectrum 
function. 

The  figures  contain  computer  produced  contours  similar 
to  those  produced  by  tne  HAIFA  code  in  previous  reports. ^ ^ 
Actual  contours  could  be  constructed  by  connecting  like  numbers? 
the  actual  value  (h  or  A)  corresponding  to  the  plotted 
numbers  are  given  by  the  formula 

Vactual  *  Vmin  +  'VVn'  [  <»♦*>«■>  1/10  . 

The  quantities  vmaxf  ^min  refer  to  maximum  arid  minimum  values 
of  the  variable  and  N  is  the  plotted  symbol.  The  ±>*  indi¬ 
cates  the  range  in  Vactual  corresponding  to  a  particular 
symbol.  The  blank  fields  appearing  in  the  plots  correspond  to 
contour  levels  intermediate  between  the  symbols  which  are  left 
blank  for  visual  clarity.  A(0,0)  appearing  on  the  spectrum 
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Figure  8  San  Joaquin-Sierra  Nevada  topography  contour  map  derived  from 
5  resolution  data. 
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Figure  9  Spectrum  function  distribution  using  topography  data 
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Figure  11  Spectrum  function  distribution  using  topography  data  from 
Figure  10. 
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function  plots  is  the  value  of  the  spectrum  function  at  k=0, 

1=0  •  Due  to  its  large  value  it  has  omitted  in  the  contouring 
so  that  resolution  could  be  improved. 

The  presence  of  the  San  Joaquin  Valley  (having  a  pro¬ 
nounced  north-west  to  south-east  axis)  in  the  grid  introduced 
an  interesting  asymmetry  to  the  topography.  Upon  examination 
of  the  graph  we  find  that  these  regularities  in  the  topography 
are  reflected  in  the  spectrum  function.  The  topography  consists 
primarily  of  two  ridges  rising  above  a  flat  plain.  One  feature 
extends  north-south  near  the  right  edge  of  the  topography 
graph  and  the  other  crosses  the  graph  at  approximately  a  45° 
angle.  These  two  ridges  correspond  to  the  two  finger-like 
extensions  seen  in  the  spectrum  function  plot  of  Figures  8 
and  10.  The  horizontal  extension  corresponds  mainly  to  the 
irregularities  in  the  N-S  ridge  while  the  diagonal  arm  reflects 
the  diagonal  ridge. 

Due  to  the  symmetry  properties  of  A,  the  graphs  exhibit 
a  reflected  symmetry  through  the  origin.  (Note:  Due  to  the 
number  of  contour  intervals  and  the  slight  displacement  of  the 
graphs  of  topography  and  spectrum  function,  exactly  the  same 
plot  character  may  not  appear  in  the  reflected  plot.)  We  also 
see  the  same  characteristic  features  appearing  in  both  A 
plots.  The  resolution  of  Figure  9  hampers  further  comparisons 
between  the  spectrum  functions  corresponding  to  5'  and  10' 
resolution.  Work  is  under  way  to  improve  this. 

Examination  of  the  other  topographic  regions  of  the 
set  reveals  that  frequently  a  prominent  topographical  feature 
will  be  reflected  in  a  clearly  observable  corresponding  feature 
of  the  A  distribution.  Further  investigation  of  these  data 
is  currently  in  progress. 
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4.5 

4.5.1 


COMPUTER  CODE  FOR  REYNOLD'S  STRESS  CALCULATION 
The  Equations 


The  equations  incorporated  in  the  code  are  substantially 
those  of  Tir  ether  ton.  [  1]  A  slight  modification  has  been  made  to 
the  trapped  wave  integrals,  however.  A  factor  of  k  in  the 
numerator  of  the  sums  which  was  omitted  by  Bretherton,  has  been 
incorporated  into  the  code. 


4.5.2  Numerical  Method 

^  Tbe  numerical  methods  discussed  in  the  previous  re¬ 
port  have  undergone  considerable  modification.  A  faster, 
more  accurate  method  of  locating  trapped  waves  was  developed 
and  incorporated.  It  was  felt  that  in  some  atmospheric  situa¬ 
tions  the  code,  in  its  present  configuration,  might  miss  a 
wave.  Another  consideration  was  the  optimization  of  the  code 
such  that  a  large  scale  survey  study  could  be  carried  out  more 
economically.  Calculations  to  date  indicate  that  a  reduction 
of  the  cost  of  a  calculation  by  a  factor  of  two  or  more  has 
been  achieved  as  a  result  of  these  steps.  Even  greater  savings 
are  expected  once  the  code  modification  program  is  completed. 

In  addition,  the  calculation  is  considerably  more  self- 
contained  and  accurate.  The  logic  of  the  new  version  of  the 
code  is  outlined  in  the  flow  chart  of  Figure  12.  Since  the  code 
is  radically  different  in  structure  from  the  previous  version, 

a  detailed  description  is  also  presented  in  the  following  sec¬ 
tion. 


A  trapezoidal  integration  scheme  is  currently  being 
employed  in  the  evaluation  of  the  stress  integrals.  The  inte¬ 
gration  uses  a  constant  increment  Ak  and  A<f>  of  k  and  <f> 

The  integration  over  k  is  currently  carried  out  for  each  value 
of  <p  ,  although  the  order  of  calculations  will  be  changed  for 
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Figure  12  General  flow  diagram  of  drag  code. 
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Figure  12  contd 
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the  newly  organized  code.  First,  the  Scorer  parameter  profile 
is  calculated  for  the  current  value  of  <J>  .  Values  of  are 
then  successively  calculated  for  all  K-values.  The  values  at 
z=0  are  then  tested  to  determine  whether  there  are  nearby 
trapped  waves.  If  the  test  is  negative  and  k  <  £,(H)  then 
the  region  is  defined  to  be  part  of  the  continuous  spectrum 
and  the  stress  integrals  are  incremented  using  Bretherton's 
Eq.  (54).  If  the  criterion  for  a  trapped  wave  is  met,  a  search 
to  localize  it  more  accurately  is  begun.  K-values  between 
k-Ak  and  k  are  searched  until  the  denominator  of  the  F 
expression  is  minimized  to  within  a  user-prescribed  accuracy. 

If  k  >  X,  (h)  ,  the  trapped  wave  occurs  at  a  K-value  for  which 
the  perturbation  velocity  w  at  z=0  vanishes.  Since  the 
real  and  imaginary  parts  of  the  solution  for  <  >  l (h)  are 
proportional  to  each  other,  they  both  must  vanish  together. 

This  property  simplifies  the  search  scheme  for  the  trapped 
wave  considerably. 

Figure  13  depicts  the  typical  K-dependence  of  the 
Scorer  equation  solutions  at  z=0  .  Let  <i  be  the  K-value 
of  the  trapped  wave.  We  also  denote  by  w1  and  w2  the  real 
parts  of  the  Scorer  equation  solution  for  k-Ak  and  k  at 
z=0  .  An  iterative  search  is  performed  which  utilizes  linear 
interpolation  between  w1  and  w2  to  obtain  a  more  accurate 
guess  at  ^  .  Using  this  value  of  k  and  its  corresponding 
wR(0)  another  linear  interpolation  is  performed  using  the 
member  of  the  old  pair  of  w' s  of  opposite  sign  tc  the  new 
value.  The  search  converges  quite  rapidly  and  usually  requires 
only  a  few  iterations.  The  criterion  for  convergence  of  the 
search  is  a  user-specified  value  of  the  permissible  change  in 
the  "s"  term  of  Bretherton's  Eq.  (61).  it  was  found  to  be 
superior  to  a  test  of  wr(k,4>)  compared  with  a  small  number 
which  was  found  to  be  unreliable  due  to  possible  rounding  errors 
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in  w  where  the  amplitude  of  w  is  very  large.  Care  must  be 
taken  in  choosing  the  size  of  Ak  since  it  might  be  possible 
for  w(k ,<p)  to  change  sign  twice  within  the  interval  and  a 
trapped  wave  would  be  missed.  Present  Ak  are  very  conserva¬ 
tive,  and  it  is  felt  that  this  situation  will  not  arise  very 
often. 

For  values  of  k  <_  £(H)  ,  the  real  and  imaginary  parts 
of  the  solution  to  Scorer's  equation  are  linearly  independent. 
Consequently,  the  real  and  imaginary  parts  will  not  simultane¬ 
ously  vanish  at  z=0  for  any  K-value.  Both  may  simultaneously 
become  quite  small,  however.  When  this  condition  occurs,  a 
peak  in  the  F  distribution  will  be  present.  The  height  of 
the  peak  is  a  function  of  w*  (0)  w(0)  ,  and  its  "Q"  ,  corre¬ 
sponding  to  the  width  of  the  peak,  is  a  function  of  the  depth 
of  the  region  in  which  the  solution  of  the  Scorer  equation 
undergoes  exponential  decrease.  A  large  value  of  "Q"  corre¬ 
sponds  to  trapped  waves  with  a  very  slight  upwards  energy  leak. 

Consequently,  a  trapped  wave  in  the  region  of  k  £  £(H) 
must  satisfy  two  conditions.  First,  the  value  of  w*(0)  w(0) 
must  have  a  local  minimum.  And  second,  the  peak  in  the  value 
of  F  must  be  sharply  localized.  Bretherton  suggests  that  the 
first  condition  can  be  recognized  by  examining  the  quantity 
arg  [w(k,  <J>)  ]  for  sudden  changes  as  a  function  of  k  .  This,  in 
effect,  requires  either  the  real  or  the  imaginary  part  of  w 
at  z=0  to  change  sign.  This  is  a  necessary  condition  for  a 
trapped  wave  to  exist.  For  <  >  Z (H )  it  is  also  a  sufficient 
condition.  However,  for  k  £  H(H)  the  second  condition  dis¬ 
cussed  above  must  also  be  met.  A  sharply  localized  peak  of  F 
requires  that  a  resonant  behavior  of  w  occur.  This,  in  turn, 
calles  for  a  sufficiently  deep  atmospheric  layer  in  which 
gravity  waves  are  reflected  so  that  the  leakage  of  energy  into 
the  stratosphere  is  small.  This  condition  can  be  quantitied  in 
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terms  of  the  quantity 


9 


which  is  a  measure  of  the  exponential  decrease  of  w  due  to 
wave  reflection.  When 


-  k2  dz  >  cr 


where  a  is  a  user-specified  number,  the  resonance  in  F  is 
judged  to  be  sufficiently  sharp  to  perform  the  integration  over 
it  analytically.  If  the  above  condition  is  not  met,  the  inte¬ 
gration  is  to  be  performed  numerically.  In  this  case,  the 
integrand  changes  sufficiently  slowly  that  an  accurate  numeri¬ 
cal  quadrature  can  be  performed  with  a  reasonable  integration 
step  Ak  . 

Since  the  extent  and  strength  of  the  region  of  exponen¬ 
tial  decay  in  the  solution  to  Scorer's  equation  increases  with 
increasing  k  ,  the  resonances  (if  any)  lying  below  k=£(H) 

exhibit  sharper  and  sharper  peaks  as  <  increases.  Through 
the  choice  of  c  the  user  has  control  over  the  extent  of  the 
spectrum  to  be  considered  as  continuum.  A  small  value  of  a 
will  require  that  rather  broad  lines  are  treated  analytically, 
while  a  large  value  will  treat  quite  sharp  lines  as  continuous. 

Tests  of  this  section  of  the  code  will  be  carried  out 
for  several  choices  of  atmospheric  parameters  in  order  to 
determine  an  optimum  value  of  a  for  accuracy  and  speed  of 
calculation. 
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5,  FUTURE  PLANS 


Now  that  the  computer  codes  for  evaluating  the  linear 
steady-state  Reynold's  stress  are  essentially  complete,  atten¬ 
tion  will  be  turned  to  the  calculation  and  parameterization 
phases.  A  full  scale  survey  study  will  involve  several  hundred 
topographical  areas  and  several  dozen  typical  atmospheric  con¬ 
ditions.  Since  each  combination  of  an  atmosphere  with  a  topo¬ 
graphic  area  constitutes  a  particular  case,  it  is  seen  that 
the  total  number  of  cases  under  consideration  is  the  product 
of  the  number  of  atmospheres  and  the  number  of  topographies. 

This  number  may  exceed  several  thousand. 

In  order  to  perform  these  calculations  in  an  efficient 
manner,  we  have  examined  the  time  requirements  of  the  codes  and 
have  considered  the  organization  of  the  calculational  sequence. 
The  linear  form  of  the  equations  is  the  key  to  the  most  ef¬ 
ficient  organization.  In  order  to  avoid  duplicate  calculation 
of  topography  factors  and  atmospheric  factors,  it  is  desirable 
to  modularize  the  steady-state  code.  Rather  than  calculating 
the  topographic  factor,  atmospheric  factor,  and  Reynold's  stress 
integral  on  a  case-by-case  basis,  it  is  desirable  to  form  and 
tabulate  the  topographic  and  atmospheric  factors  separately. 

The  desired  stress  integrals  can  then  be  found  by  combining  the 
factors  in  an  inexpensive  numerical  integration.  This  methoc* 
results  in  major  savings  in  computer  cost. 
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The  spectrum  function  for  a  typical  32x32  point  topog¬ 
raphy  requires  a  calculation  on  the  UNIVAC  1108  of  10  seconds 
CPU  time.  The  calculation  of  a  typical  atmospheric  factor  F 
requires  1100  seconds  on  the  same  computer.  For  the  survey 
study  of  approximately  100  topographies  and  10  atmospheres, 
the  case-by-case  calculation  would  require  computer  CPU  time 
of  approximately  10 6  seconds.  A  modular  approach  would  re¬ 
quire  approximately  104  seconds  for  the  atmospheres  and  10 3 
seconds  for  the  topographies.  Preliminary  estimates  are  that 
the  Reynold's  stress  integrations  will  also  require  approxi¬ 
mately  10 4  seconds.  The  total  time  would  then  be  a,  2  x  io4 
seconds,  of  which  only  a  j  nor  part  is  required  for  the  topog¬ 
raphy  Fourier  transforms. 

The  modulization  revisions  are  currently  being  made 
and  will  be  completed  shortly.  In  more  detail,  the  linear 
steady-state  calculation  is  divided  into  three  segments.  Com¬ 
munication  between  the  segments  is  implemented  through  the  use 
of  information  files  on  "fast"  drums.  The  first  segment  has 
been  completed  and  calculates  and  stores  the  spectrum  function 
distribution  for  a  given  grid.  This  code  is  described  more 
fully  in  Section  4.  The  second  segment  calculates  the  topo¬ 
graphic  factor,  F,  for  a  complete  range  of  values  of  k,<J>  . 

The  resulting  distribution  is  stored  on  mass  storage  along 
with  the  corresponding  discrete  coordinate  values.  This  code 
is  essentially  complete  and  is  substantially  the  same  as  the 
Reynold's  stress  code  described  in  Section  4.  Finally,  the 
third  segment  utilizes  the  results  of  the  first  two  to  perform 
the  integration  of  the  stress  integrals.  Since  all  of  the 
data  are  available  simultaneously,  a  higher  order  quadrature 
scheme  can  be  utilized. 

This  organization  of  the  calculation  offers  several 
advantages.  Only  specifically  desired  cases  need  be  integrated 
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and  the  files  can  be  augmented  with  additional  cases  as  de¬ 
sired.  The  individual  segments  of  the  code  can  also  be  exe¬ 
cuted  for  a  lower  computing  rate  since  a  much  smaller  portion 
of  core  is  required  than  for  the  code  version  in  which  all 
parts  are  unified.  The  computer  instructions  for  each  segment 
require  about  10,000  words  of  core  storage.  The  storage  file 
for  the  spectrum  function  requires  several  thousand  words  per 
grid.  The  F  factor  file  requires  around  10,000  words  per 
atmosphere.  From  these  numbers  it  can  be  estimated  that  the 
mass  storage  requirements  will  not  be  expensive.  Computation 
times  will  be  slightly  longer  due  to  accessing  the  drums,  but 
overall  rates  will  be  cheaper  on  a  case-by-case  comparison. 
Greater  flexibility  is  also  achieved  since  additional  atmos¬ 
pheres  or  topographies  can  be  easily  added. 
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1.  INTRODUCTION 


The  previous  project  report111  presented  the  deriva¬ 
tion  of  the  equations  solved  in  the  newly-developed  3-D 
Boussinesq  approximation  fluid  dynamics  computer  code  STUFF. 
Modeling  capabilities  were  discussed  and  to  some  extent  the 
numerical  approach  taken  in  STUFF  was  examined.  Since  then, 
the  code  has  undergone  considerable  modification.  Efforts 
have  been  primarily  directed  toward  code  optimization  in  exe¬ 
cution  in  order  to  increase  computation  speed.  However,  some 
attention  was  given  to  reducing  the  core  requirements  of  the 
instructions. 

The  following  discussion  records  the  results  of  this 
modification  program.  Also,  a  more  detailed  account  of  the 
numerical  techniques  now  employed  in  STUFF  is  given.  This 
latest  version  of  the  computer  code  is  designated  as  STUFF3. 
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2.  MODIFICATIONS  TO  THE  STUFF3  CODE 


2.1  CODE  OPTIMIZATION 

2.1.1  Introduction 

At  the  close  of  the  last  contract  period,  it  was 
determined  that  several  parts  of  the  STUFF3  code  might  be 
significantly  speeded  up  by  more  efficient  programming.  Tim¬ 
ing  runs  had  indicated  that  the  major  contributors  to  total 
execution  time  were  the  local  macro-scale  calculation  and  the 
particle  buffering  scheme.  As  a  result,  the  routines  which 
transfer  the  particles  to  and  from  core  were  completely  re¬ 
written.  The  macro-scale  calculation  was  more  efficiently 
coded,  and  a  scheme  to  selectively  avoid  the  full  macro-scale 
integration  for  some  time  steps  was  introduced.  The  results 
of  this  program  are  presented  below. 

2.1.2  Modifications  -  The  Buffering  Scheme 

To  review,  STUFF3  utilizes  a  particle  buffering  scheme 
such  that  only  a  small  fraction  of  the  Lagrangian  particle 
properties  occupy  core  storage  at  any  one  time.  The  rest  re¬ 
side  on  mass  storage  devices  in  distinct  batches.  They  are 
retrieved,  when  required,  recalculated,  and  subsequently  re¬ 
turned,  to  mass  storage.  This  method  allows  a  certain  degree 
of  machine  independence,  since  the  core  requirements  of  the 
code  are  reduced  drastically.  Currently,  STUFF 3  has  been 
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utilizing  a  particle  density  of  eight  particles  per  Eulerian 
cell.  This  results  in  a  requirement  of  64,000  particles  in  a 
20x20x20  cell  grid.  Since  each  particle  requires  six  storage 
locations  for  position,  density,  contaminant  and  turbulent 
energy  concentration,  one  sees  that  this  amounts  to  384,000 
required  storage  locations.  To  keep  such  a  large  amount  of 
information  in  core  would  be  prohibitively  expensive. 

With  the  advent  of  mass  storage  devices  with  high  rates 
of  information  transfer,  buffering  schemes  have  become  practi¬ 
cal,  both  from  an  economical  and  run  time  standpoint.  STUFF3 
stores  the  particle  variables,  associated  with  500  particles,  in 
core  at  one  time.  As  a  result  only  3000  core  locations  are 
required  to  handle  the  particle-based  information.  Since  most 
accounting  systems  utilized  by  computer  installations  heavily 
weight  core  utilization,  this  technique  results  in  marked 
economies.  The  penalty  one  pays  is  in  the  resultant  increased 
run  times  and  a  high  use  of  the  I/O  channels  to  the  peripheral 

devices.  Both  of  these  factors  increase  the  cost  of  a  calcula¬ 
tion. 

Timing  runs  were  made  which  indicated  that  the  buffer¬ 
ing  scheme  was,  indeed,  adding  substantially  to  run  costs.  As 
a  result,  an  examination  of  the  code  was  made  to  determine  if 
the  particles  contained  in  core  were  being  utilized  to  their 
maximum  efficiency.  By  alteration  of  the  order  of  the  particle 
calculation  routines,  it  was  found  that  the  number  of  exchanges 
of  particles  between  core  and  mass  storage  per  time  step  could 
be  reduced  from  9.5  to  3.5.  The  effect  of  this  reduction  has 
been  to  render  run  costs  considerably  less  dependent  on  the 
buffering  feature;  the  contribution  currently  is  approximately 
15  percent.  The  reason  for  this  small  value  is  that  arithmetic 
calculations  now  dominate  the  time  required  to  execute  a  time 
step.  Actual  run  time  and  cost  have  been  reduced  by  more  than 
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30  percent  on  an  average  basis.  Larger  grids  yield  larger 
relative  savings. 

2,1«3  Modifications  -  The  Macro-Scale  Calculation 

Further  optimization  was  centered  around  the  calcula¬ 
tion  of  the  local  turbulent  macro-scale  length.  Since  this 
requires  three-dimensional  integrations  over  all  space  and 
calculation  of  the  components  of  the  strain-rate  tensor  for 
every  cell,  a  large  fraction  of  the  total  run  time  can  be  at¬ 
tributed  to  it.  Two  optimizing  mechanisms  were  employed.  The 
first  was  a  rather  straightforward  implementation  of  a  sug¬ 
gestion  presented  in  the  last  report,  ^  namely,  that  the 
characteristics  of  the  local  mean  flow  are  influenced  primarily 
by  the  quantities  in  the  neighborhood  of  the  point  in  question, 
and  less  by  quantities  further  away.  This  scheme  was  incorpo¬ 
rated  into  the  optimization  program  as  a  limitation  on  the 
number  of  surrounding  cells  to  be  included  in  the  spatial  inte¬ 
gration.  This  limit  is  established  by  the  user,  as  an  input 
number. 

The  second  scheme  is  suggested  by  the  fact  that  the 
code  is  used  to  study  the  evolution  of  fluid  flow  from  a  non¬ 
steady  to  steady-state  configuration.  As  the  fluid  progresses 
through  the  evolutionary  sequence,  one  would  expect  that  the 
local  macro-scale  also  reach  a  steady-state  since  they  are, 
in  effect,  dependent  on  the  mean  flow.  Also,  we  would  not 
expect  large  changes  in  the  eddy  viscosity  distribution  unless 
the  macro-scales  have  changed  significantly.  The  user  speci¬ 
fies  a  value  which  he  considers  to  represent  a  significant 
change  from  one  time  step  to  the  next  in  the  macro-scale  dis¬ 
tribution.  If  this  value  is  exceeded  on  any  time  step,  then  a 
complete  re-integration  occurs.  If  not,  then  the  eddy  viscosi¬ 
ties  are  boosted  at  a  rate  which  corresponds  to  the  rate  of 
change  of  the  average  macro-scale.  This  procedure  has  the  ’ 
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advantage  of  entirely  eliminating  the  integration  at  late 

times  resulting  in  correspondingly  more  rapid  executions  of 
each  time  step. 

2-1-4  Modifications  -  Reduced  Instructions 

Secondary  optimization  focused  on  the  reduction  of  the 
number  of  machine  language  instructions  necessary  to  represent 
the  FORTRAN  coding.  To  effect  this,  all  routines  were  examined 
or  FORTRAN  structures  which  could  be  re-coded  such  that  they 
would  be  translated  by  the  compiler  more  efficiently.  This  re¬ 
coding  would  lead  to  an  overall  reduction  in  the  number  of  core 
locations  necessary  for  code  instructions.  Upon  completion  of 
this  task,  the  code  was  reduced  by  nearly  5000  words  to  18,000 
words,  a  reduction  of  over  20  percent. 

2 • 2  NUMERICAL  METHOD 
2«2.1  The  Equations 

Before  describing  the  numerical  techniques  employed  in 
STUFF3,  it  is  useful  to  summarize  the  set  of  equations  presented 

“  Prevl°us  «P°rt.  These  ar0  the  flui(J  dynaI„ic  ^ 

lence  equations  governing  the  flow.  They  are,  however,  in  a 
somewhat  awkward  form  for  numerical  solution.  A  secondary  set 
is,  therefore,  presented  which  corresponds  to  the  approach 
incorporated  into  STUFF.  The  generalized  set  is: 
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Equations  (1-3)  represent  conservation  of  mass,  mean 
flow  momentum,  and  turbulent  kinetic  energy,  respectively.  The 
fourth  equation  represents  the  thermal  energy  transport  equa¬ 
tion  written  in  terms  of  a  density  variation.  Additionally, 
with  a  redefinition  of  Q  and  D  ,  Eq.  (4)  provides  for  the 
tracing  of  an  arbitrary  passive  contaminant.  All  terms  have 
been  defined  in  the  previous  report.  The  heuristic  coeffici¬ 
ents  a,  0,  and  y  are  given  by 


)ll 
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The  local  turbulent  macro-scale  is  given  by 
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where 
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and  w(x,x')  is  a  normalized  Gaussian  weighting  function  de¬ 
fined  as 
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Additionally, 


(n1)2  =  Ofi/8xi)  Ofi/aXjL)  . 

While  Eq.  (2)  does,  in  fact,  govern  the  evolution  of  momentum 
distribution  in  the  system,  it  contains  a  pressure  gradient 
term  which  is  awkward  to  determine.  As  a  result,  many  fluid 
dynamic  computer  codes  use  a  formulation  involving  derived 
variables  which  eliminate  the  direct  pressure  calculation. 
Methods  using  the  stream  function,  for  example,  require  special 
care  in  applying  boundary  conditions,  particularly  at  internal 
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boundaries.  A  method  of  Chorin^  is  free  of  these  difficul¬ 
ties  while  still  avoiding  the  calculation  of  the  pressure  di¬ 
rectly.  The  resultant  set  is  expressed  in  terms  of  primitive 
variables,  and,  hence,  allows  straightforward  application  of 
boundary  conditions.  Equations  (1)  and  (2)  are  combined,  and 
use  is  made  of  the  assumption  that  a  pressure  field  alone  will 
not  generate  rotational  flow.  The  resulting  momentum  and  mass 
conservation  equations  appear  as: 


!t<5i>  +  3 


(V+«X/2E) 


+  Qg±  +  s 


(14) 
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With  this  modification,  the  equations  are  in  the  final  form  for 
computation.  It  now  remains  to  discuss  the  numerical  method. 

2.2.2  Numerical  Method  —  Overview 

STUFF 3  is  a  mixed-mode  hydrodynamics  code  in  which  both 
an  Eulerian  and  Lagrangian  matrix  representation  exist  concur¬ 
rently.  In  general,  mixed-mode  codes  attempt  to  minimize  the 
inadequacies  of  one  representation  by  combining  the  best  fea¬ 
tures  of  both.  This  usually  results  in  some  or  all  of  the 
field  variables  being  carried  in  both  representations  (i.e., 
available  in  either  particle  or  cell  based  quantities) ,  with 
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one  or  the  other  being  utilized  as  the  solution  progresses,  to 
evaluate  the  various  terms  appearing  in  the  differential  equa¬ 
tions.  At  some  point  there  must  exist  an  excha-ge  of  data  be¬ 
tween  the  two  such  that  each  may  be  able  to  exert  a  "correc¬ 
tive"  function  upon  the  other.  It  is  through  this  mechanism 
that  distortions  in  the  Lagrangian  net,  such  as  those  due  to 
high  amplitude  waves,  can  be  minimized  and  that  artificial 
diffusion  in  the  Eulerian  sense  is  greatly  reduced. 

Before  examining  the  data  exchange  routines  of  the 
STUFF  code,  a  discussion  of  the  allocation  of  the  dynamic  vari¬ 
ables  between  the  two  representations  is  in  order. 

All  of  the  scalar  fields  (density,  contaminant  concen¬ 
tration,  and  turbulent  energy)  are  assigned  to  both  represen¬ 
tations,  while  the  vector  fields  (in  our  case  this  is  just  the 
velocity  field)  are  associated  with  the  Eulerian  grid  alone. 

All  changes  in  the  distribution  of  a  scalar  dynamic  variable 
are  reflected  in  changes  in  the  particle  based  values  during 
a  time  step.  At  the  completion  of  the  evaluation  of  all  dif¬ 
ferential  terms,  the  Eulerian  quantities  are  assigned  the 
"average"  of  the  values  associated  with  particles  contained 
within  the  cell.  This  censusing  of  the  particles  is  one  of  two 
exchanges  of  data  in  the  STUFF  code.  The  purpose  is  primarily 
to  replace  the  explicit  Eulerian  evaluation  of  the  scalar 
advection  terms  with  a  more  accurate  Lagrangian  one.  In  this 
manner,  the  errors  induced  by  artificial  diffusion  are  avoided. 

This  censusing  procedure  is  defined  more  explicitly  in  Sec¬ 
tion  2.2.3. 

The  second  exchange  of  data  takes  place  in  the  calcula¬ 
tion  of  the  scalar  diffusion  terms.  These  terms  use  both 
particle  and  cell  based  values.  This  procedure  is  also  de¬ 
fined  more  explicitly  in  Section  2.2.3. 


SSS-R-74-2023 


2*2.3  Numerical  Method  -  A  Closer  Look 

Perhaps  the  best  method  to  describe  the  numerical  tech¬ 
niques  of  STUFF 3  is  to  detail  the  steps  of  a  typical  cycle  of 
calculation. 

Leaving  the  question  of  initialization  to  Section  2.2.4, 
we  assume  that  sufficient  data  are  available  to  begin  a  time 
step.  The  calculational  sequence  is  as  follows: 


STEP 

1 

Calculate  the  local  eddy  viscosity. 

STEP 

2 

Evaluate  source  and  diffusion  terms  in  the 
scalar  equations. 

STEP 

3 

Partially  update  particle-based  scalar  vari¬ 
ables  by  adding  source  and  diffusion  terms. 

STEP 

4 

Solve  the  momentum  equations  utilizing 
method  of  Chorin. 

STEP 

5 

Move  the  particles  with  the  new  velocity 
field  and  evaluate  the  scalar  advection 
terms . 

STEP 

6 

Remove  or  add  particles  as  needed  to  account 
for  sources  or  sinks  of  momentum. 

STEP 

7 

Establish  new  cell-based  values  from  a 
particle  census. 

STEP 

8 

Advance  the  time. 

STEP 

9 

Honor  output  requests. 

STEP  1  -  Calculation  of  the  local  eddy  viscosity. 

The  local  eddy  viscosity  is  calculated  in  routine 
EVMAKE .  Additionally,  the  turbulent  energy  source  terms  in¬ 
volving  Q  and  J  are  calculated  and  stored  in  temporary 
storage.  The  diffusion-limited  time  step  is  also  calculated. 
Finally,  a  "boosted"  eddy  viscosity  is  constructed  for  use  in 
the  momentum  equations.  This  quantity  is  the  larger  of  the 
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true  eddy  viscosity  and  a  "stability"  viscosity.  The  latter 
is  the  smallest  viscosity  needed  to  stabilize  the  weak  insta¬ 
bility  arising  from  the  finite  difference  operator  used  in  the 
momentum  equation.  The  sole  purpose  of  this  procedure  is  to 
insure  computational  stability. 

The  calculation  of  the  eddy  viscosity  requires  a  know¬ 
ledge  of  the  local  macro-scale  distribution.  Hence,  EVMAKE 
calculates  this  as  well.  At  this  point  a  feature  of  the  pre¬ 
viously  described  code  optimization  is  introduced.  A  value  is 
assigned  by  the  user  which  is  a  measure  of  the  error  that  will 
be  tolerated  in  the  eddy  viscosity  calculation.  The  average 
macro-scale  of  the  grid  is  calculated  on  each  cycle.  In  order 
for  the  local  macro-scale  to  be  spatially  integrated,  the  cur¬ 
rent  average  macro-scale  must  differ  from  the  previous  one  by  the 
given  significance  level.  If  the  test  fails,  then  the  old 
local  eddy  viscosities  are  boosted  by  an  amount  corresponding 
to  the  rate  of  change  in  the  average  macro-scale.  This  pro¬ 
cedure  is  found  to  give  excellent  results  in  practice,  showing 
good  agreement  with  the  continuous  integration  procedure.  The 
largest  deviation  appears  in  regions  of  abrupt  changes  in  the 
boundary,  as  would  be  expected. 

The  user  may  also  specify  the  maximum  extent  of  the 
spatial  integration  of  I2  and  J2  .  This  allows  substantial 
reduction  in  computing  time  with  only  a  moderate  loss  in  accu¬ 
racy  since  the  integrands  fall  off  strongly  at  larger  distances. 

STEP  2  -  Evaluation  of  scalar  source  and  diffusion  terms. 

In  this  calculation,  the  cell-based  scalar  field  values 
and  eddy  viscosity  distribution  are  utilized  to  construct  the 
diffusion  term  contributions  to  the  scalar  equations.  To  these 
terms  are  added  any  scalar  source  terms  arising  from  sources  or 
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sinks.  In  the  case  of  turbulent  kinetic  energy,  the  sources 
calculated  in  routine  EVMAKE  are  included.  The  result  is  to 
create  an  effective  "source"  rate  for  each  of  the  cells.  This 
is  stored  temporarily  for  use  in  the  next  step. 


STEP  3  -  Partial  update  of  particle  scalar  quantities. 

Some  difficulty  is  encountered  in  using  the  diffusion 
and  source  rates  to  update  the  particle  scalar  quantities  for 
a  time  At  .  The  diffusion  and  source  rates  are  cell-centered 
quantities.  A  first  approximation  would  be  to  apply  a  single 
rate  to  all  particles  falling  within  the  confines  of  one  cell. 

This  procedure  leads  to  computational  difficulties  as  seen  in 
the  figure  below. 

The  ordinate  represents  the  local  value  of  the  scalar 
variable  while  the  absissa  is  marked  off  in  cellular  intervals. 
The  dashed  lines  represent  cell-based  values  of  the  scalar 
while  the  solid  line  traces  through  the  particle-based  distri¬ 
bution.  The  arrows  represent  the  increase  or  decrease  of 
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particle-based  values  when  the  net  "source"  rates  are  applied 
to  the  particles  falling  within  each  of  the  cells.  One  sees 
that  in  many  instances  a  universal  application  of  cell-based 
rates  to  the  internal  particles  could  result  in  some  particle- 
based  values  going  negative.  As  a  result  of  this,  a  weighting 
technique  was  developed  which  drives  the  average  particle  value 
to  what  would  be  the  new  Eulerian  value  but  also  drives  the 
particles  which  deviate  most  from  the  norm  at  the  fastest  rate. 
The  effect  is  to  smooth  particle  distributions  such  that  large 
variations  are  damped  and  hence  the  Lagrangian  distribution  re¬ 
mains  undistorted  over  longer  integration  times. 

STEP  4  —  Solution  of  the  momentum  equations. 

Utilizing  the  boosted  eddy  viscosity  calculated  in 
Step  1  and  the  old  velocity  field,  the  subroutine  MOMENT 
determines  the  resultant  new  field.  The  procedure  is  to  first 
calculate  the  diffusion  terms  and  store  these  in  temporary 
storage.  Next,  the  adveotion  terms  are  calculated  and  tempo¬ 
rarily  stored.  Finally,  the  buoyancy  term  is  determined.  A 
resultant  tentative  velocity  field  is  then  determined  by 
applying  the  calculated  rates  for  a  time  At  .  This  velocity 
field  is  then  utilized  in  the  Poisson  equation  solution  for 
the  corrector  potential.  Finally,  the  resulting  corrector 
potential  is  utilized  to  update  the  tentative  velocities  to  the 
true  new  velocities  as  in  Eq.  (16).  In  order  to  achieve  a 
second  order  accurate  time  differencing  scheme  in  the  momentum 
equation  solution,  a  second  iteration  on  the  above  procedure  is 
carried  out  with  the  "old"  velocity  field  being  replaced  by  an 
average  of  the  true  old  velocity  field  and  the  one  just  calcu¬ 
lated.  it  is  found  that  this  second  iteration  is  not  very 
costly  in  terms  of  increasing  total  cycle  execution  times. 

This  is  not  obvious  since  the  Poisson  equation  solution  requires 
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an  iterative  procedure  in  its  solution.  It  turns  out  that  an 
over-relaxed  Gauss-Seidel  procedure  produces  quite  rapid  con¬ 
vergence.  Additionally,  a  conservative  value  of  the  Courant 
number  also  aids  in  rapid  convergence. 

STEP  5  —  Move  the  particles. 

Once  the  new  velocity  field  is  determined,  the  evalu¬ 
ation  of  the  scalar  advection  terms  is  carried  out.  This  is 
done  by  moving  the  particles  at  their  local  velocity  for  a 
time  At  .  The  local  velocity  is  determined  from  two  consid¬ 
erations:  first,  the  ime  centered  velocity  field  is  used, 
and  second,  an  interpolation  scheme  is  used  to  determine  the 
velocity. appropriate  to  the  position  of  the  particle. 

STEP  6  —  Add  or  delete  particles  to  account  for  sources  of 
momentum. —  - - 

Lagrangian  marker  particles  must  be  added  or  removed 
as  necessary  from  the  calculation  as  the  fluid  element  in 
which  they  reside  is  swept  into  or  out  of  the  computational 
mesh.  Additionally,  the  discrete  representation  of  the  veloc¬ 
ity  field  will  sometimes  result  in  particles  being  moved  into 
internal  obstacle  cells  accidently.  These  conditions  are 
handled  in  two  ways  depending  on  the  boundary  conditions. 

If  the  cyclic  boundary  option  is  specified,  then 
Par^^c^-es  swept  out  of  the  downstream  side  are  reinserted  at 
the  upstream  side  and  complete  their  movement  for  that  time 
step.  Any  particle  accidentally  being  moved  into  an  obstacle 
region  is  placed  at  the  position  it  occupied  before  the  move¬ 
ment.  This  insures  a  constant  number  of  particles  within  the 
fluid. 

If  an  a  priori  boundary  condition  is  specified,  then 
Par^^c^-es  which  are  swept  out  of  the  downstream  side  are 
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flagged  as  deleted  in  their  respective  batches.  The  vacant 
space  then  becomes  available  for  the  insertion  of  new  particles 
generated  by  fluid  entering  on  the  upstream  side,  or  from  in¬ 
terior  sources  of  momentum.  New  particles  are  given  scalar 
values  characteristic  of  the  source  from  which  they  emerged. 
Presently,  the  code  deletes  all  particles  which  are  acciden¬ 
tally  moved  into  an  interior  obstacle  when  the  a  priori  bound¬ 
ary  condition  is  specified.  This  tends  to  reduce  the  number 
of  active  particles  slightly  but  produces  little  alteration 
of  results  when  compared  to  the  repositioning  method. 

All  of  these  calculations  occur  in  routine  PMOVE  which 
also  performs  the  search  of  the  particle  batches  for  particles 
flagged  as  deleted.  It  then  packs  the  batches  eliminating  the 
deleted  particles  and  inserting  any  newly  generated  ones  into 
the  vacancies. 

STEP  7  —  Establish  new  cell-based  values. 

With  the  completion  of  the  movement  of  the  particles, 
the  equations  effectively  have  been  solved  at  the  new  time. 

The  task  remains,  however,  to  generate  the  new  cell-based 
quantities.  These  quantities  are  required  in  the  calculation 
of  additional  time  steps  and  are  desired  in  the  editing  of 
problem  results. 

The  cell-based  quantities  are  calculated  in  routine 
CENSUS  which  utilizes  essentially  the  same  procedure  found  in 
the  MAC  method.  Each  particle  is  visualized  to  extend 
Ax^/N  in  each  of  the  cellular  dimensions  and  to  have  a  scalar 
value  throughout  this  volume  equal  to  the  value  attributed  to 
the  particle.  N  is  the  cube  root  of  the  particle  density. 

The  result  is  to  smooth  the  distribution  and  to  reduce  the 
likelihood  of  large  discontinuities  arising  from  particle 
clumping.  The  CENSUS  routine  determines  the  fraction  each 
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P  rticle  donates  to  each  cell.  These  contributions  are  summed 
and  normalized  to  yield  a  resultant  cell-based  distribution. 
Finally,  CENSUS  checks  to  insure  that  no  cell  is  empty  of 
particles,  if  any  are  found  in  this  condition,  then  the  cell 
value  is  given  the  averaged  value  of  all  the  adjacent  cells. 
This  rarely  happens  in  runs  where  the  particle  density  is 
eight  per  cell  or  greater. 

STEP  8  —  Advance  the  time. 

The  time  is  now  advanced  so  that  a  new  cycle  can  begin. 
Several  stability  criteria  representing  different  terms  in  L 
equations  limit  the  size  of  the  time  step. 

in  order  to  avoid  the  problems  associated  with  having 
o  introduce  more  than  one  set  of  source  particles  per  time 
step  across  any  boundary  of  the  domain,  the  stability  number 
is  restricted  to  values  of  f  <  1/N  .  This  allows  consider¬ 
able  coding  economy  and  does  not  restrict  the  allowable  time 
steps  too  severely. 

The  stability  criterion  in  explicit  schemes  requires 
a  restriction  on  At  due  to  diffusion.  This  is  calculated  in 
routine  EVMAKE  and  has  the  following  form: 


At  < 


47  l 


-+-L-  + 

Ax2  Ay2 


1 

Az2 


where  Ax,  Ay,  Az  represent  cellular  dimensions,  e  is  the 
local  eddy  viscosity,  and  f  is  a  stability  number  such  that 

,,  An°ther  stability  criterion  incorporated  into  STUFF 
leots  the  stability  associated  with  the  propagation  of  in- 

r!quireWaVeS'  ^  ^  6nSUre  Stabilit>'  in  this  case,  we  • 
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At  <  fAx/c 

where  Ax  is  the  smallest  dimension  of  the  cell,  f  a  stabil¬ 
ity  number,  and  c  the  internal  wave  phase  velocity  given  by 

where  X  is  the  largest  wavelength  of  the  generated  internal 
waves,  g  the  gravitational  acceleration  and  A£  the  maximum 
fractional  variation  in  density  across  the  cell  interfaces. 

X  ,  in  effect,  is  a  free  parameter  to  be  chosen  a  priori  from 
past  experience  with  waves  arising  in  fluid  flow  problems. 

Its  value  is  not  critical  for  most  cases. 

Finally,  the  stability  condition  for  the  advection 
scheme  can  be  expressed  as  follows: 


un+1  + 


.n+1 . 2 


At  < 


r  +  2f  Ax 


Account  has  been  taken  above  of  the  possible  increase  in  veloc¬ 
ity  during  the  anticipated  time  step. 

The  final  resultant  time  step  to  be  used  in  the  next 
increment  is  chosen  to  be  the  most  restrictive  of  the  four  in¬ 
equalities  presented  above. 


STEP  9  —  Honor  output  requests. 

Requests  for  edits  of  output  data  can  assume  three 
forms:  (1)  requests  for  a  printer  dump  of  all  variables  of 

interest  inclu^i^a  vector  and  scalar  arrays,  (2)  requests  for 
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printer  plots  of  the  velocity  and  scalar  fields;  an  option  is 
to  be  specified  by  the  user  of  which  particular  variables  are 
desired.  The  plots  are  depicted  as  if  one  were  looking  down 
from  the  top  onto  the  model  and  successively  removing  layers  of 
cells.  Finally,  (3)  a  tape  dump  may  be  made  for  the  purposes 
of  saving  the  variable^  for  later  use,  such  as  a  restart. 

Other  output  features  were  mentioned  in  the  previous  report. ^ 

This  completes  a  calculational  cycle  in  STUFF3.  We 
next  discuss  the  initialization  of  a  problem.  The  options 
available  for  starting  a  calculation  are  outlined  below. 

2.2.4  Initialization 

The  Velocity  Field  —  Three  methods  are  available  to 
the  user  in  specifying  an  initial  configuration.  The  first, 
which  is  of  limited  use,  involves  the  determination  of  the  in¬ 
ternal  momentum  distribution  from  potential  flow  and  specified 
boundary  flow.  This  method  requires  the  assumption  that  the 
boundary  flow  must  be  irrotational.  This  corresponds  to  an 
input  flow  uniform  in  z  with  no  differential  cross  flow,  a 
condition  not  likely  to  be  found  in  atmospheric  studies. 

A  second  technique  is  to  specify  the  velocity  in  the 
x "direction  to  be  periodic  and  proceed  as  above.  This  has  the 
advantage  of  uniquely  determining  the  flow.  A  proper  boundary 
condition  is  also  imposed  at  the  downstream  grid  face  which 
does  not  introduce  flow  disturbances.  However,  this  technique 
is  of  doubtful  usefulness  in  3-D  codes  since  large  numbers  of 
downwind  cells  are  required  so  that  wrap  around  effects  do  not 
become  significant.  This  approach  becomes  prohibitively  ex¬ 
pensive  in  realistic  situations. 

The  third  method  appears  to  be  the  most  useful  for 
mountain  wave  problems.  This  involves  the  useage  of  the  cor¬ 
rector  potential-velocity  relations  expressed  by  Eqs.  (15)  and 
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(16) .  The  procedure  is  to  use  the  best  guess  of  the  ambient 
field  available.  This  will  be  in  considerable  error  near 
highly  convoluted  regions  and  the  divergence  will  probably  be 
non- zero.  This  field  is  then  utilized  as  a  source  term  of 
the  Poisson  equation  to  derive  a  corrector  potential  which, 
when  combined  as  in  Eq.  (15)  with  the  guessed  ambient  field, 
will  result  in  a  field  which  is  non-divergent  and  preserves 
the  vorticity  distribution  of  the  original  guess.  This  pro¬ 
cedure  gives  realistic  ambient  fields,  and  although  they  are 
not  unique,  they  have  represented  excellent  starting  fields 
for  several  test  cases.  This  method  will  allow  the  specifi¬ 
cation  of  essentially  arbitrary  z-dependent  u— velocity  pro¬ 
files.  Cross  winds  ere  also  allowed. 

The  Scalar  Fields  — —  Initialization  of  the  Eulerian— 

based  values  of  the  scalar  fields  was  described  in  the  pre- 

f  11 

ceedmg  semiannual  report.  1  The  particle-based  values  are 
initialized  from  the  above  Eulerian-defined  distributions. 
Essentially,  every  particle  which  falls  within  a  given  cell 
is  initialized  with  the  values  of  the  Eulerian-based  quantities 
of  density,  contaminant  concentration,  and  turbulent  energy 
density. 
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3.  CONCLUSIONS 

With  the  description  of  initialization  procedures  we 
have  essentially  completed  our  discussion  of  the  Boussinesq 
hydro— code  STUFF3.  Some  additional  detail  regarding  particular 
topics  is  presented  in  the  appendices. 

STUFF3  represents  one  of  the  more  advanced  fluid  dy¬ 
namics  codes  in  use  in  atmospheric  research  study.  The  mixed¬ 
mode  structure  and  incorporation  of  a  heuristic  model  of  tur¬ 
bulence,  both  represent  an  excellent  utilization  of  current 
state-of-the-art  knowledge.  Machine  requirements  are  not 
significantly  greater  than  previous  Eulerian  codes  and  execu¬ 
tion  times  are  comparable. 

Future  use  of  the  code  is  to  be  discussed  in  the 
"Future  Plans"  section. 
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4.  FUTURE  PLANS 


The  above  discussion  essentially  completes  the  code 
development  work  on  the  3-D  transient  model.  Its  future  ap¬ 
plication  will  be  limited  due  to  higher  priority  commitments, 
especially  the  wave  drag  parameterization  study.  If  time  and 
funds  permit,  a  limited  3-D  transient  wave  drag  study  using 
STUFF3  will  be  conducted. 

The  2-D  transient  model  HAIFA  will  be  used  as  needed 
in  the  parameterization  study  to  resolve  questions  regarding 
transient  effects.  STUFF 3  will  also  be  utilized  for  this  pur 
pose. 
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APPENDIX  A 


THE  GRID 


The  grid  utilized  by  STUFF3  is  pictured  below.  For 
computational  efficiency  and  boundary  condition  application, 
the  true  grid  is  surrounded  on  all  sides  by  a  layer  of  border 
cells.  The  velocities  at  the  interior  faces  of  these  cells 
can  be  assigned  to  apply  various  boundary  conditions  on  the 
flow.  Work  is  almost  complete  to  incorporate  outflow  boundary 
conditions  utilizing  extrapolated  internal  values.  This  pro¬ 
cedure  is  equivalent  to  the  setting  of  the  second  derivative 
of  the  normal  component  to  zero,  thereby  minimizing  the  dis¬ 
turbing  effect  of  the  boundary  on  the  internal  flow. 

Boundary  cells  and  interior  obstacle  cells  are  flagged 
with  a  value  of  zero  in  the  grid-defining  array,  FULL(I,J,K)  . 
Any  cell  characterized  by  zero  (meaning  zero  fluid)  is  assigned 
velocities  associated  with  its  interfaces  at  problem  generation 
time.  This  allows  interior  obstacle  cells  to  act  as  ducts 
through  which  fluid  can  be  removed  from  or  injected  into  the 
system,  i.e.,  these  become  sources  or  sinks  of  momentum. 

The  particles  are  initially  distributed  within  the  grid 
at  a  specified  density  (in  terms  of  number  of  particles  (N)  per 
cellular  dimension).  Hence,  a  three-dimensional  problem  will 
have  N3  particles  per  cell.  The  initial  spatial  distribution 
within  the  fluid  filled  cells  is  seen  in  Figure  A-2. 
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APPENDIX  B 

DEFINING  GRID-BASED  VARIABLES 

The  indexing  notation  used  to  access  the  various  scalar 
and  vector  quantities  is  shown  for  an  arbitrary  i,j,k  cell  in 
Figure  F-l.  All  scalar  quantities  are  cell-centered  with 
i,j,k  r aferring  to  the  cell  center  coordinates.  Velocities 
are  interface  quantities  with  the  u,v,w  components  of  a  cell 
being  defined  on  the  i+^j,  j+*j,  and  k+*s  faces,  respectively. 
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Figure  B-l  Notation  utilized  by  STUFF.  An  xy  cross-section 
is  shown  taken  through  the  cell  center.  Scalars 
are  cell-centered  quantities.  Vector  components 
are  interface  quantities. 
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APPENDIX  C 

PROCEDURE  IN  UPDATING  PARTICLE  ARRAYS 
DUE  TO  SOURCE  AND  DIFFUSION  TERMS 


In  a  purely  Eulerian  scheme  the  updating  of  the  field 
variables  due  to  diffusive  contributions  is  straightforward. 
Utilizing  the  standard  centered  difference  operator,  the  cell- 
based  quantity  tends  with  time  to  a  value  more  in  equilibrium 
with  its  neighbors.  Since  STUFF  utilizes  the  Eulerian-based 
diffusion  rates  to  modify  the  particle-based  values  in  the 
cell,  it  is  possible,  using  an  average  rate,  to  force  values 
associated  with  some  of  the  particles,  negative.  While  a  sub¬ 
sequent  census  would  yield  a  cell-based  value  in  agreement 
with  the  purely  Eulerian  method,  discontinuities  may  develop 
when  these  particles  undergo  subsequent  advection.  For  ex¬ 
ample,  some  of  the  negative  particles  might  be  moved  into 
cells  in  which  t.iey  are  the  sole  contributors.  A  subsequent 
census  would  then  yield  negative  cell-based  quantities.  In 
order  to  avoid  this  problem  a  smootning  function  was  developed 
which  preferentially  modifies  particle  values  and  reduces  the 
variation  from  the  expected  Eulerian  value. 

This  function  should  be  related  to  the  local  curvature 
in  the  field.  Also,  it  should  be  functionally  dependent  on 
the  strength  of  the  diffusive  coefficient  and  the  time  interval 
over  which  it  can  act.  After  much  experimentation,  the  follow¬ 
ing  form  was  selected: 
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0P±1  =  nn  .  /nn  _  n  \  s-AtS» 

WPART  yCELL  \UPART  QCELL/ 


Where  Q  represents  the  field  variable,  PART  and  CELL  refer 
to  particle  and  cell-based  values,  respective! y .  E  represents 
the  sum  of  the  turbulent  eddy  viscosity  and  molecular  diffusion 
Coefficients,  and  S  is  related  to  the  field  uhape  by: 


Where  the  summation  range  over  all  cells,  and  the  derivatives 
fire  approximated  by  standard  cell  centered  finite  difference 
operators.  The  last  term  involving  repronents  the  rate  of 

phange  due  to  diffusion  and  sources  and  is  dot  ormined  from  all 
based  quantities  at  t  =  tn  as  is  the  S  factor. 
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APPENDIX  D 

THE  PARTICLE  CENSUSING  SCHEME 


Once  the  particles  have  been  moved  to  their  final 
positions  for  the  current  time  step,  a  census  is  carried  out 
to  determine  the  resultant  cell-based  values  so  that  a  new 
time  step  may  commence.  For  the  purposes  of  the  census,  ench 
particle  is  considered  to  fill  a  finite  volume,  whose  dimen¬ 
sions  are  given  by  Ax±  =  AX^N  .  is  the  particle  dimen¬ 

sion  in  the  ith  direction,  AX±  is  the  cell-dimension  in  Uie 
ith  direction,  and  N  is  the  cube  root  of  the  particle  den¬ 
sity  per  cell.  This  description  is  defined  graphically  in 
Figure  D-l. 

The  calculation  of  the  cell-based  quantities  now  im- 
comes  a  series  of  summations,  whereby  the  value  for  a  particu¬ 
lar  cell  is  given  by: 


<P 


PELL 


dv. 

l 


where  $  is  the  scalar  variable  with 
ring  to  particle  and  cell-based  values, 
is  the  part  of  the  volume  of  particle  i 


PART  and  CELL  refcr- 
respectively.  dv^ 
falling  within  the 


boundaries  of  the  given  cell. 
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Figure  D 


1  A  cross-section  of  Eulerian  cells  containing 
a  particle  "cube".  In  this  example,  particle 
density  is  8/cell,  hence  Ax.  =  Ax./2  . 
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APPENDIX  E 

THE  PARTICLE  MOVEMENT  SCHEME 


The  scalar  advection  terms  are  calculated  through  the 
actual  moving  of  the  particles.  After  the  completion  of  this 
procedure  a  census  is  taken  to  establish  the  new  cell-based 
values.  This  results  in  an  evaluation  of  the  advection  terms 
which  is  free  of  artificial  diffusion. 

Some  errors  are  introduced  by  this  scheme,  however. 
While  they  are  much  less  significant  than  those  found  in  the 
method  they  replace,  they  should  be  noted. 

Each  particle  in  the  fluid  should  be  advected  by  the 
fluid.  In  order  for  this  to  be  the  case,  the  instantaneous 
velocity  of  the  particle  must  be  the  local  velocity  of  the 
fluid.  A  finite  difference  scheme  can  only  approximate  the 
velocity  field  at  each  space-time  point.  As  a  result,  the 
velocity  which  they  experience  is  some  interpolated  value 
between  locally  defined  values.  Once  this  value  has  been 
determined  the  particle  is  moved  at  this  velocity  for  the 
time  At  . 

Two  types  of  errors  are  introduced.  The  first  results 
from  the  fact  that  the  calculated  instantaneous  velocity  of 
the  particle  is  an  interpolated  one.  The  resulting  error  will 
be  small  if  the  Eulerian  grid  is  of  the  proper  resolution  to 
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adequately  resolve  the  problem.  A  second  error  is  introduced 
when  the  particle  is  moved  at  the  constant  velocity  for  a  time 
At  .  Actually,  the  local  velocity  is  a  function  of  t  .  The 
magnitude  of  this  error  is  still  significantly  less  than  the 
scheme  being  replaced,  however. 

In  order  to  insure  that  the  scheme  is  second  order 
accurate  in  time,  the  time  centered  velocity  field  is  utilized 
for  the  particle  movement. 


